%Conduct a variety of spectral analyses, evaluate the goodness of fit
%for H0, and give the indicated statistical significance at the
%orbital bands.

clear;
addpath ../Tools;  
%Spectral toolbox described in Huybers, Putting the significance of spectral
%peaks on the level: implications for the 1470-yr peak in Greenland
%18O, Journal of Climate, 2022.  

Params.nw=[3 4 5]; 
Params.p=[1 2];
Params.q=[0 1 2];
Params.dt=[1 5]; 
Params.diff=[0 1];
Params.qredo=1; 

if Params.qredo,
    load ../Bathymetry/collect_all.mat All;
    count=0;
    for ct1=1:length(Params.nw),
        for ct2=1:length(Params.p),
            for ct3=1:length(Params.q),
                for ct4=1:length(Params.dt),
                    for ct5=1:length(Params.diff),
                        count=count+1;
                        Pall(count).setting=[ct1 ct2 ct3 ct4 ct5];
                        [count ct1 ct2 ct3 ct4 ct5],
                        opts.nw=Params.nw(ct1);
                        opts.p=Params.p(ct2);
                        opts.q=Params.q(ct3);
                        opts.dt=Params.dt(ct4);
                        opts.diff=Params.diff(ct5);
                        
                        dof=2*(2*Params.nw(ct1)-1);
                        Vexp(count)=psi(1,dof/2); 

                        for ct=1:length(All.S);
                            if rem(ct,25)==0, ct, end;
                            z=All.z(ct,:);
                            z(isnan(z))=nanmean(z);
                            z=smoothPH(z,opts.dt);
                            z=z(1:opts.dt:end);
                            if opts.diff==1, z=diff(z); end;
                            [P s coeffs]=specw(detrend(z),opts);
                            Pall(count).Vexp=Vexp(count); 
                            Pall(count).P(ct,:)=P;
                            Pall(count).s=s;
                            Pall(count).V(ct)=sum(P.^2)/length(P); 
                        end;
                        return;
                    end;
                end;
            end;
        end;
    end;
    save full_spectral_analysis.mat;
else
    load full_spectral_analysis.mat;
end;
%Results of this analysis are shown in Table S2.  

for ct=1:length(Pall); 
    count(ct)=ct;
    Vexp(ct)=Pall(ct).Vexp;
    nw(ct)  =Params.nw(Pall(ct).setting(1));
    p(ct)   =Params.p(Pall(ct).setting(2));
    q(ct)   =Params.q(Pall(ct).setting(3));
    dt(ct)  =Params.dt(Pall(ct).setting(4));
    diff(ct)=Params.diff(Pall(ct).setting(5));

    %significance of orbital bands
    clear P;
    for ct2=1:length(Pall(ct).P(:,1));
        P(ct2,:)=exp(Pall(ct).P(ct2,:));
    end;
    
    pl=find(All.S>4);
    k=length(pl)*2*nw(ct)-1;
    Pm41=nanmean(P(pl,:)); 
    [dummy plf]=min(abs(Pall(ct).s-1/41));
    pval41(ct)=1-gamcdf(Pm41(plf),k,1/k);
    Vavg(ct)=nanmean(Pall(ct).V(pl)); 
end;

X=[count', nw', p', q', dt', diff', pval41', Vexp', Vavg'];
matrix2latex(X,'pvals.tex','format', '%-6.3f');

%Summary statistics:
sum(All.S>4),

%number significant
[length(pval41), sum(pval41<0.01), sum(pval41<0.05)],

%number significant that are within 20% of expected F
F=Vavg./Vexp;
pl=find(F>=0.8 & F<=1.2);
[length(pl), sum(pval41(pl)<0.01), sum(pval41(pl)<0.05)],

%formulation insignificant
[dummy pl]=max(pval41);
[pl F(pl) pval41(pl) nw(pl) p(pl) q(pl) diff(pl)],

%formulation shown in Figure 3
pl=find(diff==0);
[dummy pl2]=min(F(pl));
pl=pl(pl2); 
[pl F(pl) pval41(pl) nw(pl) p(pl) q(pl) diff(pl)],

%formulation with minimun F
[minF pl]=min(F);
[minF nw(pl) p(pl) q(pl) diff(pl)],

%formulation with maximum F
[maxF pl]=max(F);
[maxF nw(pl) p(pl) q(pl) diff(pl)],